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An infinite projected entangled-pair state (iPEPS) is a variational tensor network ansatz for 2D 
wave functions in the thermodynamic limit where the accuracy can be systematically controlled 
by the bond dimension D. We show that for the doped Hubbard model in the strongly correlated 
regime {U/t = 8, n = 0.875) iPEPS yields lower variational energies than state-of-the-art variational 
methods in the large 2D limit, which demonstrates the competitiveness of the method. In order 
to obtain an accurate estimate of the energy in the exact infinite D limit we introduce and test 
an extrapolation technique based on a truncation error computed in the iPEPS imaginary time 
evolution algorithm. The extrapolated energies are compared with accurate quantum Monte Carlo 
results at half filling and with various other methods in the doped, strongly correlated regime. 
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I. INTRODUCTION 

The accurate study of strongly correlated systems is 
one of the biggest challenges in condensed matter physics. 
A well known example is the 2D Hubbard modet^ which 
potentially captures the relevant physics of the cuprate 
high-Tc superconductors. Despite its simplicity and an 
enormous effort in trying to solve the model, the phase di¬ 
agram of the Hubbard model is still controversial. Still, 
in recent years substantial progress has been achieved 
with a variety of different numerical methods (see e.g. 
Refs. HQ, so that there is hope that the full solution 
of the Hubbard model may become within reach in the 
near future. For recent state-of-the-art benchmark re¬ 
sults from various methods, see Ref. H 

Solving systems in ID is far better under control than 
in 2D, mostly thanks to the well-known density-matrix 
renormalization group (DMRG) method.^i^ DMRG has 
an underlying variational tensor network ansatz, called 
matrix product state (MRS), in which the wave function 
is efficiently represented by a trace over product of ten¬ 
sors. The accuracy of the ansatz can be systematically 
controlled by the bond dimension D of the tensors, and 
one typically reaches extremely accurate results in ID 
(and quasi ID). DMRG can also be used to study 2D 
systems (typically on cylinders) by mapping th e sy stem 
onto a ID problem with long-ranged interactions.^^ How¬ 
ever, the computational cost scales exponentially with 
the width of the cylinder such that the approach is not 
scalable to large 2D systems.^ 

In order to overcome this exponential scaling 2D ten¬ 
sor network ansatze have been developed, such as pro¬ 
jected entangled- pair st ates (PEPS,'i^l^ also called ten¬ 
sor product state^i^R^ or th e 2D multi-scale entangle¬ 
ment renormalization ansatz.These networks are 
designed in such a way that they reproduce an area 
law scaling of the entanglement entropy which a large 
class of relevant ground states in 2D fulfill.^ The in¬ 


volved methods are technically more complicated than 
MPS-based approaches which is one of the main rea¬ 
sons why it took several years to develop these methods. 
However, recently there have been substantial break¬ 
throughs which clearly demonstrate the enormous poten¬ 
tial of 2D tensor networks. For example, it was shown 
for the t- J modeP^ that infinite PEPS (iPEPS) - an 
ansatz for a state in the thermodynamic limit - yields 
lower variational energies than the state-of-the-art re¬ 
sults from fixed-node Monte Carlo.l^ Another example 
is the Shastry-Sutherland model in a magnetic field,!^ 
where iPEPS helped to gain a new understanding of the 
magnetization process, thanks to largely unbiased simu¬ 
lations. 

Thus, already current (i)PEPS algorithms can outper¬ 
form (or compete with) the best variational methods for 
strongly correlated fermionic models like the t-J model 
or also for frustrated spin systems (see e.g. Refs. [25l - 
dZI. In this paper we show that the same is true also in 
the strongly correlated regime of the 2D Hubbard model, 
where we find lower variational energies than the best 
variational Monte Garlo results for large 2D systems.!^ 

One major difficulty in iPEPS simulations so far has 
been to obtain an accurate estimate of the energy in the 
exact, infinite D limit. Typically the energy does not 
smoothly depend on the bond dimension D, making an 
extrapolation of the finite D data to the infinite D limit 
problematic. Accurate extrapolations become particu¬ 
larly important if several states at finite D are strongly 
competing, as e.g. in the t-J model,^ where uniform and 
stripe states exhibit almost the same energy at finite D. 
A precise estimate of the energy is crucial to identify the 
true ground state among these competing states. 

In this paper we propose and test an approach to ex¬ 
trapolate the energy based on a truncation error w which 
quantifies the degree of approximation in the iPEPS 
imaginary time evolution algorithm. This quantity plays 
a similar role as the truncation error in conventional 
DMRG simulations which is typically used to extrapo- 
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late energies. Empirically we find here that the energy 
varies in a much smoother way with w than with 1/D 
such that an extrapolation in re —>■ 0 yields an improved 
estimate of the exact ground state energy. 

We benchmark this extrapolation technique for the 2D 
Hubbard model, first in the exactly solvable U/t = 0 
limit (which is particularly challenging for iPEPS since 
the ground state is strongly entangled), then at finite U/t 
at half-filling where we compare our results with accurate 
Quantum Monte Carlo results.l^ Finally, we also provide 
an estimate of the energy in the doped, strongly corre¬ 
lated regime {U/t = 8, n = 0.875) and a comparison with 
various other methods from Ref. [S] 

This paper is organized as follows. In the next sec¬ 
tion we give a short introduction to the iPEPS ansatz 
and the ground state algorithm based on imaginary time 
evolution. In Sec. Ill we discuss ways to perform energy 
extrapolations with iPEPS, in particular, we explain how 
to compute a truncation error and use this quantity as an 
extrapolation parameter. In Sec. IV we present our finite 
D and extrapolated energies for the 2D Hubbard model 
and a comparison with other methods. Finally, in Sec. V 
we summarize our findings and give prospects for solving 
the Hubbard model. In appendix A we present addi¬ 
tional results for the 2D Heisenberg model and a model 
of non-interacting spinful fermions with a pairing poten¬ 
tial, to provide further evidence for the usefulness of the 
extrapolation technique based on the truncation error. 


II. IPEPS ANSATZ AND METHOD 

An infinit e projected entangled-pair state 
(iPEPS)P^E^E^ is an efficient variational tensor network 
ansatz for two-dimensional states in the thermodynamic 
limit which obey an area law of the entanglement 
entropy (typically ground states of local Hamiltonians). 
The ansatz consists of a supercell of tensors which is 
periodically repeated on a lattice, with one tensor per 
lattice site.^ On the square lattice, each tensor has a 
physical index and four auxiliary indices which connect 
to the nearest-neighboring tensors. The accuracy of the 
ansatz can be systematically controlled by the bond 
dimension D of the auxiliary indices (i.e. each tensor 
contains variational parameters where d is the local 
dimension of a lattice site). An iPEPS with D — 1 
corresponds to a product state, and by increasing D 
entanglement can be added in a systematic way. 

For translational invariant states a supercell with only 
one single tensor can be used. If the translational sym¬ 
metry is spontaneously broken, a supercell compatible 
with the symmetry breaking pattern is needed (e.g. for 
an antiferromagnetic state two different tensors for the 
two sublattices are required). Since in practice the struc¬ 
ture of the ground state is not known in advance, we 
run simulations using different supercell sizes to check, 
which supercell yields the lowest variational energy. This 
approach also provides a way to determine several com¬ 


peting low-energy states, as for example done in the t-J 
model, ^ in which uniform and different stripe states have 
been found using different supercell sizes. In order to find 
the true ground state among these competing states, one 
needs to have an accurate estimate of the energy of each 
state in the infinite D limit. However, a simple extrapo¬ 
lation in D often fails to give an accurate estimate, due 
to the non-smooth dependence of the energy on D, and 
this is why it is important to find alternative ways to 
perform such extrapolations. Such an improved extrap¬ 
olation technique will be presented in the next section. 

The iPEPS wave function is evaluated by contracting 
the two-dimensional tensor network in a controlled ap- 
proximate way. In the present work we use a varia nl)^ 
of the corner-transfer matrix (CTM) method.l^SIini The 
accuracy of the contraction is controlled by the ’’bound¬ 
ary” dimension %, which we choose large enough (up to 
several hundreds) such that the resulting error is negligi¬ 
ble (compared to the effect of the finite D). To inc rease 
the efficiency we make use of abelian symmetries! ^^ * ^^ ! 
For an introduction to iPEPS, see e.g. Refs. [33] and |34l 
We note that 2D tensor networks have first been intro¬ 
duced for spin systems, and later extended to fermionic 
systems, see Refs. |33l |35Hl2| 

In order to obtain an approximate representation of 
the ground state for a given Hamiltonian iJ, the tensors 
need to be optimized, i.e. the best variational parame¬ 
ters have to be found. For iPEPS this is typically done 
by performing an imaginary time evolution of an initial 
(e.g. random) iPEPS. The evolution operator is split 
into a product of two-site operators via a Trotter-Suzuki 
decomposition (assuming nearest-neighbor interactions). 


exp(-/3i7) « Db j , Df, = exp(-Ti7b), (1) 

where the product goes over all nearest-neighbor bonds b 
in the supercell, Hf, is the Hamiltonian term on bond 5, 
and T = /3/n is a small imaginary time step. The imag¬ 
inary time evolution is performed by sequentially multi¬ 
plying the two-site operators U^ to the iPEPS and rep¬ 
resenting the resulting wave function again as an iPEPS, 
until convergence is reached. 

Let us consider the application of such a two-site op¬ 
erator Ub to two tensors A and B which are connected 
by the bond b. The resulting state \^'aib') = Ub\'^ ab) 
can be represented by two new tensors A' and B' where 
the bond dimension on bond b has increased from D to 
D' < d?D. For an efficient evolution the corresponding 
bond needs to be truncated back to the original bond di¬ 
mension D, resulting in a truncated wave function Ilk j^). 
In the so-called full updat^^ (or fast full updat^^ this 
truncation is done by finding the new tensors A and B 
which minimize the cost function 

C = = min Y^d(i,R), (2) 

A,B A,B ’ 
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with 

d{A,B) = {^'a>b'Wa'B') + {^abI^ab) (3) 
-{'^'a'B'I^ab) - {^abI'^'a'B')- 
Finding the new tensors can be solved in an iterative way, 
as explained e.g. in Refs. [551 and IMP^ 

III. ENERGY EXTRAPOLATION WITH IPEPS 

Typically, for challenging problems, one does not reach 
convergence as a function of D and one needs to perform 
an extrapolation to the infinite D limit to obtain an esti¬ 
mate of the true ground state energy. One possibility is 
to plot the energy as a function of 1/ZI and then trying 
to extrapolate the data to l/D -a 0. However, in prac¬ 
tice the energy does not depend on 1 /D in a smooth way 
which makes an accurate extrapolation difficult (see e.g. 
Refs, mi mi mi and 051 for examples). 

Empirically one finds that the overall convergence of 
the energy goes faster than linear in 1/D, such that a 
linear extrapolation in 1 /D (using the largest few values 
of D) provides a lower bound Ei of the true ground state 
energy. Since the method is variational, the energy for 
the largest value of D corresponds to an upper bound E^. 
In practice a crude estimate of the energy can be obtained 
from the mean value Em = {EuEE/)/2 with an error bar 
of A = {Eu — Ei)/2. Examples of these estimates will 
be shown in the next section. While this approach can 
provide a reasonable guess of the exact energy, a more 
accurate and controlled extrapolation would be highly 
desirable. 

In DMRG simulations energy extrapolations are typi¬ 
cally much more accurate by extrapolating in the trun¬ 
cation error e, corresponding to the sum of the dis¬ 
carded squared singular values in the two-site variational 
optimization.!^ In simple words, the truncation error 
measures how far away the state is from the true ground 
state, which is reached if e goes to zero. The question is 
now if a similar quantity could also be used in iPEPS sim¬ 
ulations to improve energy extrapolations. The most nat¬ 
ural way would be to implement a similar two-site vari¬ 
ational optimization algorithm in iPEPS. However, the 
imaginary time evolution algorithm is more commonly 
used and more easy to implement, and we therefore aim 
to define a similar quantity within this approach.!^ 

Let us consider the cost function C in Eq. ([^ . For the 
true ground state, which is reached for sufficiently large 
D and /3, the cost function C is zero. However, if the 
true ground state is not reached because D is too small, 
then the cost function will reach a certain non-zero value 
C[D,fi -A oo) for large /?. For small r the cost func¬ 
tion depends linearly on the time step (to lowest order). 
We now define the quantity w{D) = C{D,I3 -A oo)/t 
which is independent of r (for small t) and decreases 
monotonously with D. This quantity w measures the 
truncation error when approximating an iPEPS with en¬ 
larged bond dimension D' on a bond, \'^'a'B')^ with a new 


iPEPS with smaller bond dimension D, . Thus, w 

plays a similar role as the truncation error e in DMRG 
simulations. 

A priori we do not know how the energy depends on w, 
in contrast to DMRG simulations where the energy con¬ 
verges linearly in the truncation error for sufficiently large 
bond dimensions (after a suitable number of finite-system 
sweeps). Nevertheless, we find here that the energy de¬ 
pends on w in a much more regular way than on 1/D, 
such that an extrapolation in w using a polynomial fit 
to the data provides an improved estimate of the exact 
energy in the re —>■ 0 limit. This will be illustrated with 
several examples for the Hubbard model in the next sec¬ 
tion. 


IV. BENCHMARKS: 2D HUBBARD MODEL 

As a benchmark we consider the single-band 2D Hub¬ 
bard model with only nearest-neighbor hoppings, 

H = -t (cLg<t + H.c.) +UY (4) 

where cl^ (ci^) creates (annihilates) an electron with spin 
= {ti)-} on site i, fiia = is the number operator, 

U is the on-site repulsion, and t the hopping amplitude. 

In the following we first present the iPEPS results in 
the non-interacting case D/t = 0 at half filling (n = 1) 
which can be exactly solved. For the interacting case 
U/t > 0 we make a comparison with several other 
methods from the recent state-of-the-art benchmark pa¬ 
per by LeBlanc et al.,!^ including auxiliary-field quan¬ 
tum Monte Carlo (AFQMC), density-matrix embedding 
theory (DMET), DMRG, and fixed-node Monte Carlo 
(FNMC). We first consider the half filled case for {7/t = 4 
and U/t = 8 which can be accurately solved by AFQMC 
simulations since there is no sign problem at half fill¬ 
ing. Finally, results for the doped case n = 0.875 in the 
strongly correlated regime U/t = 8 are presented and 
compared to the best available data. All energies are 
given in units of t. 


A. U/t = 0 at half filling 

It is known that 2D free fermionic systems with a ID 
Fermi surface have a multiplicative logarithmic correc¬ 
tion to the area law of the entanglement entropy.!^ Since 
an iPEPS can only reproduce an area law this case poses 
a particular challenge and we do not expect to obtain 
the exact result for finite D (on an infinite lattice). Nev¬ 
ertheless, since the area law is only weakly violated (in 
contrast to a volume law), one can still obtain an approx¬ 
imation to the ground state and a variational estimate of 
the ground state energy. Eor example, for D = 16 iPEPS 
yields an energy per site oi E = —1.597. Compared to 
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FIG. 1. (Color online) iPEPS energies as a function of the in¬ 
verse bond dimension (squares) and as a function of the trun¬ 
cation error w (circles) for the non-interacting case Ult — 0. 


the exact result —1.6211..., this corresponds to a relative 
error of « 1.5%. 

The iPEPS energies as a function of 1/D are shown 
in Fig. (squares). One can clearly see how the varia¬ 
tional energy improves upon increasing D. However, the 
dependence on 1/D is not very regular which makes an 
extrapolation in l/H somewhat difficult. As discussed in 
the previous section we can obtain a lower bound, Ei, of 
the ground state energy by a linear extrapolation of the 
data in 1/D, whereas the value at the largest D corre¬ 
sponds to a variational upper bound, The range of 
energies estimated in this way is [—1.597, —1.636], which 
includes the exact value. Based on this data we obtain 
an estimate Em = —1.616 ± 0.019. 

We next consider the iPEPS data plotted as a function 
of the average truncation error w (circles in Fig. Q. One 
can clearly observe a smoother dependence in E(w) than 
in the E{l/D) data. Fitting the data with a third-order 
polynomial yields an energy Au, = —1.6217 in the limit 
w —f 0 which is close to the exact energy. As an estimate 
of the error we take half the difference between lowest 
variational energy and the extrapolated value, A = {E^ — 
E^)/2 = 0.012, shown by the dark blue error bar in 
Fig.0 As an alternative estimate we average over several 
fits using different ranges of data points, which yields 
Etu = —1.6219 with a standard deviation of 0.006, shown 
by the light blue error bar in Fig. 

Thus, even in the ’’worst-case” for iPEPS (i.e. free 
fermionic systems) we can obtain quite an accurate esti¬ 
mate of the ground state energy based on an extrapola¬ 
tion in the truncation error w. 



0 0.05 0.1 0.15 

1/D or w 


FIG. 2. (Golor online) iPEPS results for the energy for 
U/t = 4 at half filling (n = 1) as a function ol 1/D (squares) 
and w (circles), in comparison with the results from other 
methods (extrapolated to the thermodynamic limit). The 
reference value obtained from AFQMG is shown by the green 
line, with an error bar indicated by the dashed lines. 


B. U/t = 4 at half filling 

Next we consider the case U/t = 4 at half filling, where 
the ground state is insulating and exhibits antiferromag¬ 
netic long-range order. Since there is no ID Fermi surface 
as in the U/t — 0 case, we expect the ground state to be 
less entangled (obeying an area law). This is reflected 
in a higher accuracy of the ground state energy obtained 
with iPEPS, and in a smaller truncation error w com¬ 
pared to the U/t = 0 case. For example, taking D = 16 
the relative error (compared to the extrapolated AFQMC 
resuHP —0.8603 ± 0.0002) is of the order of 0.14%, i.e. 
an order of magnitude better than the U/t = 0 result. 
The truncation error is w{D = 16) ^ 0.05 compared to 
w{D = 16) ^ 0.16 in the U/t = 0 case. 

The iPEPS energies exhibit a rather irregular behavior 
as a function of 1/D, shown by the squares in Fig. A 
linear fit using the five largest D values yields the lower 
bound El = —0.8610. Computing an estimate based 
on the 1/D extrapolation as in the U/t = 0 case yields 
—0.8602 ± 0.0008 in agreement with the AFQMC result. 

The energies as a function of w (circles) exhibit a more 
regular behavior, as previously observed in the U/t = 
0 case. A third order polynomial fit including all data 
points yields Eyj = —0.8603±0.0005. If we average again 
over several fits using different ranges of data points we 
obtain E^ = —0.8604 ± 0.0005 in agreement with the 
AFQMC result. 

Our data is also in agreement with the extrapolated 
DMET and DMRG results, shown in Fig. with a com¬ 
parable error bar. 
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FIG. 3. (Color online) iPEPS results for the energy for 
I//t = 8 at half filling (n = 1) in comparison with other meth¬ 
ods. 


C. U/t = 8 at half filling 

As we move away from the non-interacting limit, 
iPEPS becomes more accurate (in contrast to weak- 
coupling approaches). For U/t = 8, D = 16 the energy is 
—0.52415. The relative error compared to the AFQMC 
result (—0.5247±0.0002) is small, only 0.1%. Thus, even 
without using any extrapolation we obtain already a re¬ 
markably accurate result in the thermodynamic limit. 

The iPEPS estimate in the infinite D limit based on the 
1/D extrapolation is —0.5246±0.0005 (using the 4 largest 
D values) and —0.5250 ± 0.0008 (using the 5 largest D 
values), in agreement with the AFQMC result. 

Also here, the E{w) curve is much smoother than the 
E{l/D) data, shown in Fig. A third-order polyno¬ 
mial fit yields Eyj = —0.5244 ± 0.0001, which is slightly 
higher, but still compatible with the AFQMC result. (A 
similar result is obtained by taking the average over sev¬ 
eral fits.) 

Compared to the other methods iPEPS shows the best 
agreement with the AFQMC data: DMRC is slightly too 
high (—0.5241 ± 0.0001), DMET has a large error bar 
(—0.5234 ± 0.001) and the FNMC estimate is too high 
(-0.52315 ± 0.00005). 


D. Doped case for U/t = S 

Finally we present results in the doped, strongly cor¬ 
related regime for U/t = 8 and n = 0.875 in Fig. |^for 
which AFQMC is no longer exact due to a strong sign 
problem. Determining the ground state in this regime 
is one of the fundamental open problems in the field 
of strongly correlated systems. The main reason why 


this is so challenging is that there are many different 
states which are energetically very strongly competing 
(e.g. uniform superconducting states or different types 
of stripe states), which explains why different state-of- 
the-art methods yield different answers for the physics in 
this regime, see e.g. Ref. [9l 

One of the main advantages of iPEPS is that the differ¬ 
ent competing states can be obtained by using different 
supercell sizes, and the properties of these states (e.g. 
order parameters) can then be studied individually for 
each low-energy state. This was done for example in 
Ref. \n\ for the t-J model (an effective model of the Hub¬ 
bard model in the strongly interacting limit) where e.g. 
a uniform d-wave superconducting state with coexisting 
antiferromagnetic order is found by using a 2-site super¬ 
cell, or period-5 stripe in a 5 x 2 supercell. In order to 
identify the true ground state among these competing 
states, one needs to compare their extrapolated energies, 
and this why an accurate extrapolation of the energy is 
important. 

The systematic study of all possible supercell sizes and 
the discussion of the physics of all the competing states 
is beyond the scope of this paper and left for future work. 
Instead we illustrate here the advantage of the improved 
extrapolation technique by comparing the energies of two 
different competing states: A vertical stripe state with 
a period 5 with charge and spin orders coexisting with 
superconductivity, and a (non-superconducting) diagonal 
stripe state with a period 16 with one hole per unit length 
per stripe. 

If we consider the energies plotted as a function oll/D, 
shown in Fig. we can observe that for a given D the 
diagonal stripe state has a lower variational energy than 
the vertical stripe state. However, since the energies of 
both states are still rapidly decreasing with increasing D 
it is hard to predict which state is lower in energy in the 
large D limit based on a 1/D extrapolation, which yields 
Em = -0.763 ± 0.010 and E^ = -0.764 ± 0.015 for the 
vertical and diagonal stripe, respectively. These values 
lie very close and have a large overlapping error bar. 

If, however, we use the truncation error w as an ex¬ 
trapolation parameter, we can obtain a much clearer dis¬ 
tinction between the two states. The extrapolation in 
w yields E^ = -0.7637 ± 0.005 {E^ = -0.7577 ± 0.004) 
for the vertical (diagonal) stripe; averaging over several 
fits using different ranges yields = —0.7633 ± 0.002 
[Ew = —0.7581 ± 0.0014). Thus, in the large D limit the 
vertical stripe state is favored over the diagonal stripe 
state. 

We next compare our results for the vertical stripe 
states with other methods from Ref. m Our best iPEPS 
variational energy for D = 16 is E = —0.75325. This 
is lower than the best variational result —0.74884 from 
FNMC for a 20 X 20 system where the FNMC energies 
are increasing with system size. Thus, iPEPS clearly pro¬ 
vides a lower variational energy in the thermodynamic 
limit than FNMC, which demonstrates the competitive¬ 
ness of iPEPS in the doped, strongly correlated regime. 
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FIG. 4. (Color online) iPEPS energy of a vertical stripe with 
period 5 and a diagonal stripe with period 16 for the doped 
Hubbard model in the strongly correlated regime {U/t — 8, 
n = 0.875) in comparison with other methods. 


The extrapolated iPEPS energy is comparable to 
the result by (approximate) constrained path AFQMC 
(—0.766 ± 0.001), the DMRG result (extrapolated in 
the truncation error) for a finite cylinder of width 6 
(—0.759 ± 0.004), and the DMET result in a 5 x 2 cell 
(-0.7671).Eal 

V. SUMMARY AND PROSPECTS 

We presented iPEPS results for the energy of the 2D 
Hubbard model at half filling for U/t = 0, 4, and 8, and 
away from half filling for U/t = 8, n = 0.875. In the 
doped case the variational energies at large bond dimen¬ 
sion are lower than the best available variational Monte 
Carlo results which demonstrates that iPEPS is a very 
competitive variational method in the strongly correlated 
regime. It is in this challenging and physically relevant 
region where iPEPS (or 2D tensor networks in general) 
have the largest potential to go substantially beyond the 
present state-of-the-art. 

In order to obtain an estimate of the exact ground state 
energies we proposed to perform an extrapolation in the 
truncation error w, complementary to the (more crude) 
extrapolations in I/D, allowing us to compute ground 
state energies with a higher accuracy. At half filling the 
extrapolated results agree with the exact value in the 
non-interacting case, and with accurate AFQMC results 
for {7/t = 4 and U/t = 8. These extrapolations will play 
a key role to identify the true ground state among sev¬ 
eral competing states (e.g. stripe and uniform state^^ 
which lie very close in energy. As an example we pro¬ 
vided an estimate of the energy of a (vertical) period-5 
stripe for U/t = 8, n = 0.875, and showed that the ex¬ 


trapolated energy is lower than the one of a period-16 
diagonal stripe. 

The present iPEPS data has been obtained using mod¬ 
est computational resources. We believe that with large 
scale parallel simulations using bond dimensions up to 
D ~ 20 ... 24, in combination with the present extrapo¬ 
lation technique, the ground state phase diagram in the 
strongly correlated regime {U/t > 6) is accessible. Com¬ 
bined with appro ach es which work best in the weakly 
correlated regime,^^ and with supporting results from 
other methods in the strongly correlated regime, the full 
solution of the 2D Hubbard model seems within reach. 
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Appendix A: Additional results 

In order to further demonstrate the usefulness of the 
energy extrapolation technique based on the truncation 
error w we present additional results for the 2D Heisen¬ 
berg model and for an exactly solvable model of non¬ 
interacting spinful fermions with a pairing potential in 
this appendix. 


1. 2D Heisenberg model 

We consider the two-dimensional S=l/2 Heisenberg 
model on a square lattice with Hamiltonian, 

H = J^SiSj, (AI) 

(ul) 

where Si is a spin-1/2 operator on site i. We set the cou¬ 
pling J = 1. Since there is no negative sign problem this 
model can be solved by Quantum Monte Carlo, and ac¬ 
curate estimates of the energy can by obtained by an ex¬ 
trapolation of the finite-size data to the thermodynamic 
limit. The Monte Carlo estimate is E = —0.669437(5) 
from Ref. [50] obtained from linear system sizes up to 
L — 16. A more precise estimate was presented in 
Ref. [Sn E = —0.6694421(4), using larger system sizes. 

Also in this example the iPEPS energies as a function 
of 1/D show a rather irregular behavior, as shown in 
Fig. Using an 1/D extrapolation yields —0.66943(6) 
(using the four largest values of D), in agreement with 
the extrapolated QMC result. 
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FIG. 5. (Color online) iPEPS energy per site of the 2D Heisen¬ 
berg model compared to the Quantum Monte Carlo result. 
Note that the truncation error on the x-axis has been rescaled 
by a factor 5 for better visibility. 


The truncation errors in this example are considerably 
smaller than for the Hubbard model, indicating that the 
ground state is less entangled. Also in this case the de¬ 
pendence on w is much smoother than on 1/D. A sec¬ 
ond order polynomial ht including all data points yields 
Ew = —0.66945(4). A similar result is obtained by aver¬ 
aging over several fits including different ranges of data 
points. 

2. Non-interacting spinful fermions with a pairing 
potential 

In this section we consider an exactly solvable model 
of spinful fermions with a nearest-neighbor hopping and 
a pairing potential term, 

H =-t {c\^cj„ +H.c)j 

(hi.O’) 

{i,j) 


with 7 ij the amplitude of the pairing potential. In the 
present example we set t = 1, and 'jij = ±1 for bonds 
oriented along the x- and y-direction, respectively. The 
exact energy in the thermodynamic limit is —1.35494.... 

The iPEPS energies shown in Fig. [^decrease rapidly 
with increasing D. As a consequence, the linear extrapo¬ 
lation in 1/D leads to an estimate with a very large error 
range, —1.36 ± 0.01. 

In contrast, performing an extrapolation in the trunca¬ 
tion error yields a much better estimate, = —1.3547± 
0.002, or Em = —1.3547 ± 0.001 when taking the average 
over several fits. 

This example illustrates that estimating the correct 
energy based on a l/E extrapolation can be hard, even 
though the exact value lies not far from the energy at the 
largest bond dimension. Thanks to the much smoother 
behavior of the E{w) curve one can obtain a much more 
accurate estimate using an extrapolation in w. 



FIG. 6. (Color online) iPEPS energy per site of an exactly 
solvable model of spinful fermions compared to the exact re¬ 
sult. 
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